import pandas as pd
import numpy as np
import pickle
import warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
Wykorzystanym zbirem danych jest Home Equity (HMEQ), zawierający informacje o 5960 klientach banku, którzy otrzymali kredyty hipoteczne.
Na podstawie zbioru próbowałam przewidzieć prawdopodobieństwo defaultu, czyli faktu, że klient będzie zalegał z płatnościami – określa to binarna zmienna BAD (1 oznacza default). Pozostałe 12 zmiennych opisuje m.in. historię kredytową aplikującego, historię zawodową oraz charakterystyki obecnej pożyczki.
Więcej informacji na temat danych można znaleźć pod linkiem https://www.kaggle.com/ajay1735/hmeq-data
hmeq = pd.read_csv("hmeq.csv", error_bad_lines=False)
hmeq_info = {'BAD' : 'client defaulted on loan 0 = loan repaid',
"LOAN" : "Amount of the loan request",
"MORTDUE" : "Amount due on existing mortgage",
"VALUE": "Value of current property",
"REASON": "DebtCon debt consolidation HomeImp = home improvement",
"JOBS" : "occupational categories",
"YOJ": "Years at present job",
"DEROG" : "Number of major derogatory reports",
"DELINQ": "Number of delinquent credit lines",
"CLAGE": "Age of oldest trade line in months",
"NINQ": "Number of recent credit lines",
"CLNO": "Number of credit lines",
"DEBTINC" : "Debt-to-income ratio"}
from pandas.api.types import is_numeric_dtype
{column : is_numeric_dtype(hmeq[column]) for column in hmeq.columns}
set(hmeq['REASON'])
set(hmeq['JOB'])
hmeq = pd.concat([hmeq, pd.get_dummies(hmeq['REASON'], prefix='REASON', dummy_na=True)],axis=1)
hmeq = pd.concat([hmeq, pd.get_dummies(hmeq['JOB'], prefix='JOB', dummy_na=True)],axis=1)
hmeq.drop(['REASON', 'JOB'],axis=1, inplace=True)
hmeq.isna().sum()
hmeq_nonan = hmeq.dropna()
X = hmeq_nonan.iloc[:, 1:]
y = hmeq_nonan.loc[:, "BAD"]
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.35, random_state=42)
X_val, X_test, y_val, y_test = train_test_split(X_test, y_test, test_size=0.6, random_state=42)
for data in [X_train, X_test, X_val, y_train, y_val, y_test]:
data.reset_index(drop=True, inplace = True)
X_train.shape
metrics = ["accuracy_train", "accuracy_test", "roc_auc_train", "roc_auc_test"]
from sklearn.ensemble import RandomForestClassifier
rf_final1 = pickle.load(open("final_nonan_rf.p", "rb"))
from sklearn.metrics import accuracy_score, roc_auc_score
results = {metric : {} for metric in ["accuracy_test", "roc_auc_test"]}
results["accuracy_test"]["RandomForest"] = (accuracy_score(y_test, rf_final1.predict(X_test)))
results["roc_auc_test"]["RandomForest"] = (roc_auc_score(y_test, rf_final1.predict_proba(X_test)[:,1]))
results = pd.DataFrame(results)
results
np.random.seed(42)
ind = np.random.randint(len(X_train.index)+1, size=1)[0]
obs = pd.DataFrame(X_train.iloc[ind, :]).T
y_train[obs.index].values
obs.size
rf_final1.predict_proba(obs)
Przykładowe wyjaśnienie metodą breakdown.
import dalex
from dalex.explainer import Explainer
exp = Explainer(model = rf_final1, data = X_train, y = rf_final1.predict_proba(X_train)[:,1], model_type = "classification")
exp_test = Explainer(model = rf_final1, data = X_test, y = rf_final1.predict_proba(X_test)[:,1], model_type = "classification")
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot()
Przykładowe wyjaśnienie metodą SHAP.
import shap
shap.initjs()
X,y = X_train.reset_index(drop=True), y_train.reset_index(drop=True)
model = rf_final1
# explain the model's predictions using SHAP
# (same syntax works for LightGBM, CatBoost, scikit-learn and spark models)
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X, check_additivity=False)
# visualize the first prediction's explanation (use matplotlib=True to avoid Javascript)
shap.force_plot(explainer.expected_value[0], shap_values[0][0,:], X.iloc[0,:])
shap.force_plot(explainer.expected_value[0], shap_values[0], X)
important_variables = {}
for i in range(len(X_train.index)):
obs = pd.DataFrame(X_train.iloc[i, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
df = break_down.result.iloc[break_down.result["contribution"].abs().argsort()]
var_sign = df.loc[:, ["variable_name", "sign"]]
indexNames = var_sign[var_sign['variable_name'].isin(["intercept", "prediction", ""]) ].index
var_sign.drop(indexNames , inplace=True)
important_variables[i] = {"variable_name" : var_sign.loc[:, "variable_name"].values,
"sign" : var_sign.loc[:, "sign"].values}
variables = [value['variable_name'] for key, value in important_variables.items()]
signs = [value['sign'] for key, value in important_variables.items()]
variables_2 = [i[-2:] for i in variables]
signs_2 = [i[-2:] for i in signs]
Obserwacjami o różnych najbardziej wpływowych zmiennych okazały się m.in. te o indeksach 2 i 10 ze zbioru treningowego.
ind_1 = 2
ind_2 = 10
Co ciekawe, mają one te same wartości zmiennej zależnej.
print(y_train.reset_index(drop=True)[ind_1])
print(y_train.reset_index(drop=True)[ind_1])
W przypadku obserwacji o indeksie 2 są to zmienne DEBTINC i DELINQ oznaczające odpowiednio stosunek długu do przychodu i liczbę istniejących kredytów, w których aplikujący zalega ze spłatami, natomiast na wynik obserwacji 10 najbardziej wpłynęły CLAGE i MORTDUE, czyli "wiek" najstarszej linii kredytowej oraz kwota do spłacenia na istniejącym kredycie.
print(variables_2[ind_1])
print(variables_2[ind_2])
for feature in list(variables_2[ind_1])+list(variables_2[ind_2]):
print(feature + ": " + hmeq_info[feature])
Możemy przyjrzeć się tym obserwacjom na wykresach. Poniżej pierwsza z nich, czyli obserwacja o indeksie 2.
obs = pd.DataFrame(X_train.iloc[ind_1, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))
Możemy zauważyć, że wartość stosunku długu do przychodu ok. 28 zmniejsza predykcję (czyli zmniejsza prawdopodobieństwo wejścia w default) o 0.055, za to liczba istniejących kredytów równa 1 zwiększa ją o 0.11.
Poniżej wykres rozbicia obserwacji o indeksie 10.
obs = pd.DataFrame(X_train.iloc[ind_2, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))
Widzimy, że najbardziej wpływową zmienną MORTDUE, czyli kwota pozostała do spłacenia na istniejącym kredycie, która będąc równa ok. 50000 zmniejsza predykcję o 0.031. Odwrotny efekt ma zmienna CLAGE (wiek najstarszej linii kredytowej), która przy wartości 93 zwiększa predykcję o 0.026.
vars_signs = [[(k, l) for k, l in zip(i, j)] for i, j in zip(variables, signs)]
Wybraną przez mnie zmienną do porównania jest DEBTINC czyli stosunek kwoty długu do przychodu.
Obserwacjami, które będę porównywać są obserwacje ze zbioru treningowego o indeksach 0 i 1.
ind_1=0
ind_2=1
Mają one inne wartości zmiennej zależnej.
print(y_train[0])
print(y_train[1])
print(vars_signs[ind_1][-1])
print(vars_signs[ind_2][-1])
Wizualizacja rozbicia pierwszej obserwacji.
obs = pd.DataFrame(X_train.iloc[ind_1, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))
Widzimy, że wartość zmiennej DEBTINC równa 34 powoduje obniżenie predykcji o 0.022.
Wizualizacja rozbicia drugiej obserwacji.
obs = pd.DataFrame(X_train.iloc[ind_2, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))
Tutaj z kolei możemy zaobserwować, że wartość tej zmiennej równa 46 wpływa na znaczne podwyższenie predykcji.
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(solver="liblinear", random_state=42, intercept_scaling= 1.3).fit(X_train, y_train)
exp_lr = Explainer(model = lr, data = X_train, y = lr.predict_proba(X_train)[:,1], model_type = "classification")
important_variables_lr = {}
for i in range(len(X_train.index)):
obs = pd.DataFrame(X_train.iloc[i, :]).T
break_down = exp_lr.predict_parts(obs, type = "break_down")
df = break_down.result.iloc[break_down.result["contribution"].abs().argsort()]
var_sign = df.loc[:, ["variable_name", "sign"]]
indexNames = var_sign[var_sign['variable_name'].isin(["intercept", "prediction", ""]) ].index
var_sign.drop(indexNames , inplace=True)
important_variables_lr[i] = {"variable_name" : var_sign.loc[:, "variable_name"].values,
"sign" : var_sign.loc[:, "sign"].values}
variables_lr = [value['variable_name'] for key, value in important_variables_lr.items()]
signs_lr = [value['sign'] for key, value in important_variables_lr.items()]
variables_2_lr = [i[-2:] for i in variables_lr]
signs_2_lr = [i[-2:] for i in signs_lr]
ind = 0
obs = pd.DataFrame(X_train.iloc[ind, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))
obs = pd.DataFrame(X_train.iloc[ind, :]).T
break_down = exp_lr.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))
W przypadku lasów losowych najbardziej wpływową zmienną jest znany nam już DEBTINC, czyli stosunek długu do dochodu, który znacząco obniża predykcję.
W regresji logistycznej ma on dużo mniejszy wpływ (co więcej – w drugą stronę), za to na prowadzenie wychodzi zmienna YOJ oznaczająca liczbę lat w obecnej pracy, która w modelu lasów losowych miała bardzo niewielki wpływ, którego zwrot także różnił się od zwrotu w modelu regresji.
Jak zostało pokazane, wyjaśnienia poszczególnych obserwacji w obrębie modelu mogą być bardzo różne, nawet gdy wartości zmiennej zależnej są takie same. Można też było zauważyć, że choć dana wartość zmiennej obniża predykcję, niewielkie, zdawałoby się, zwiększenie tej wartości skutkuje całkowicie różnym wpływem na predykcję.
Dodatkowo, różnice mogą występować także w obrębie jednej obserwacji -- w zależności od modelu zmienne wpływające na jej predykcję mogą się bardzo różnić.